Phase transition in inelastic disks 



Teruhisa S. Komatsu 

Department of Pure and Applied Sciences, University of Tokyo, Komaha, Meguro-ku, Tokyo 153, 

Japan 

Abstract 

This letter investigates the molecular dynamics of inelastic disks without 

external forcing. By introducing a new observation frame with a rescaled 

time, we observe the virtual steady states converted from asymptotic energy 

dissipation processes. System behavior in the thermodynamic limit is carefully 

investigated. It is found that a phase transition with symmetry breaking 

occurs when the magnitude of dissipation is greater than a critical value. 
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In energy conservative systems, it is well known that macroscopic properties at an equi- 
librium state can be described by a few state variables |1|]. This raise a question as to 
whether, in dissipative systems, similar state variables might be defined for macroscopic 
properties at steady states maintained by an external energy source. When a dissipation 
rate is sufficiently small, the system might be expected to behave as a conservative system; 
but is this actually the case? To examine this and other related questions, it is interesting 
to investigate models that can connect a conservative and a dissipative system by varying 
the system parameters. 

An ensemble of elastic hard disks is one of the simplest models for describing the fiuid 
state in conservative systems. Computer simulations played an important role in discovering 
the existence of fiuid-solid transition in this model . We will here consider an ensemble of 
inelastic hard disks, which has previously been investigated as a model of granular materials 
0. By varying the inelasticity, we can continuously change the system from a conservative 
to a dissipative system. This ability will be useful for our present purposes of investigating 
a thermodynamic properties in dissipative systems. 

In order to attain a steady state in a dissipative system, an external energy source that 
compensates for energy dissipation due to inelastic collisions is indispensable. In granular 
physics, the vibrating bed is a typical energy source. Nevertheless, attaching such an 
energy source will break the isotropy of the system from the outset. In this letter, we 
simulate the dissipative system without energy sources and investigate the virtual steady 
states using a new observation frame that will be introduced later. 

The inelastic disk system without energy input , known as a cooling gas or dissipative gas 
system, has been investigated using molecular dynamics (MD) and hydrodynamic models. 
These studies reveal that the homogeneous state is unstable for sufficiently high inelasticity 
even without attaching any energy source. In this state, a cluster of particles ^ and collec- 
tive mean flow (shear mode) are formed. We study how the instability is characterized 
by means of MD simulation of an ensemble of inelastic hard disks, with a particular focus 
on the asymptotic collective phenomena and its system size dependence. 
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The system consists of an ensemble of inelastic hard disks in two dimensional space. 
For simplicity, each particle is assumed to have a unit mass and its rotational degrees of 
freedom are ignored. Collisions between circular particles are inelastic; the inelastic collision 
is implemented in the following manner. At a collision of two particles, the tangential 
velocities to the collision plane are preserved, while the normal component of the relative 
velocity Af„ changes to Af^, where 

A< = -e„At;„. (1) 

Here, the coefficient of restitution e„(0 < e„ < 1) is constant for all collisions. In the 
following, we parameterize the inelasticity by e(= 1 — e„). In the case e = 0, the system 
becomes conservative with the elastic collisions; in the case e > 0, energy dissipation occurs 
due to the inelastic collisions. 

The selected boundary conditions consist of a square box enclosed by elastic rigid walls. 
Thus the energy of the system is not dissipated by the bouncing of particles against the 
wall. For initial conditions, we adopt spatially homogeneous states equilibrated by setting 
e = 0. The time evolution of the system is calculated by the event-driven method [0]. 

In this letter, we focus on asymptotic behavior of the system after a sufficiently long 
transient time. We note that our model is almost identical with the model in |^ except for 
boundary conditions. 

Under our boundary conditions, there are three important system parameters: the num- 
ber of particles A^, the inelasticity e and the area fraction of particles 0. Here 0(0 < < 1) 
is the ratio of the total area covered by particles to that of the box S. 

In the conservative case e = 0, it is known that this system, consisting of an ensemble of 
elastic hard disks, has two phases: a fluid and a solid phase. In this letter, we investigate 
the system in the range of parameter corresponding to the fluid phase when e = 0. For 
this parameter range, our results do not depend on the dispersity of particle radii a. In the 
following, we show only the results for monodisperse case(a = 0.5); however, the results of 
our simulations for the polydisperse case (uniform distribution in the range 0.4 < a < 0.5) 
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are almost the same. On the other hand, in the range of parameter corresponding to the 
sohd phase when e = 0, the results depend on the dispersity of particle radii. The results 
for this parameter range will be reported in a future publication. 

For sufficiently large values of e, it is known that the time development by the event 
driven method is ill-defined because infinite collisions occur in a finite time interval (inelastic 
collapse) p. Our investigation is limited to sufficiently small values of e in order to avoid 
the inelastic collapse. 

We will now describe how to observe the system. Since no energy source is attached to 
the system, the total energy of the system monotonically decreases in time. Nevertheless, 
by the rescaling of time, it is possible to make the energy of the system virtually conserved. 
We introduce the rescaled time t and the rescaled velocity of i-th particle Vi. 

i = v^^*^' = Mt)/\fm, (2) 

where T(t) is the averaged kinetic energy per one particle at a time t, T{t) = J2ivf(t) /2N] 
and Vi{t) is the original velocity of i-th particle. We call the rescaled system "R-system" 
and the original system "0-system" . All of the present observations are carried out for the 
"R-system". Under the translation by Eqs.(^, the averaged kinetic energy per one particle 
in the "R-system" is normalized to the unity, T{t) = Y^i'^fit)/'^-^ — 1' i-^-' total kinetic 
energy is conserved for any e in the "R-system". Further, asymptotic energy dissipation 
processes in the "0-system" are translated to steady states in the "R-system" . 

It must be noted here that the translation by Eqs.(0) is just the rescaling of time and that 
the trajectories of the particles in space are not infiuenced by this rescaling. We calculate the 
time evolution of the "0-system" by the event driven method, and observe the "R-system" 
with a special focus on its steady states after a sufficiently long transient time. In the rest 
of this letter, all variables shown refer to those in the "R-system" . 

The dependencies of pressure on for several values of e are shown in Fig.|I|. Pressure P 
is defined as the time averaged sum of the impulses at the bouncing of particles on the walls 
per unit length per unit time in the "R-system". The vertical axis in the figure is NT /PS, 
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which is unity for the ideal gas hmit (0 — > 0). Here 5* is the area of the box. 

In the conservative case e = 0, only the fluid phase exists in the system for cj) < (pc- 
The value of (pc is reported to be (pc — 0.7 Consider the dissipative cases e > in 

Fig.|l](a) paying attention to the difference from the case e = 0. When e = 0.02, there is no 
remarkable difference. For e = 0.04, a decrease in pressure is observed in the intermediate 
range of 0. For e = 0.08, a similar difference appears in the wider range of 0. 

Comparing Fig.0(a) and (b), it is found that similar dependencies of pressure on are 
observed for the same values of A^e. Thus A^e is an important parameter of the system and 
it would characterize the distance from the conservative system. 

For sufficiently large A^e, the behavior of the system is quite different from the case e = 0. 
A snapshot of the system for (A^, e) = (256, 0.08) is shown in Fig.0(b) compared to that for 
e = (Fig.|(a)). 

From Fig.0(b), it is found that a correlation of particle velocity exists and that a mean 
flow circulating anti-clockwise in the box is formed. As a result of the emergence of the 
circulating mean flow, the impulse of particles to the walls at bounce decreases compared 
with the case e = 0. Then the decrease of pressure, i.e., the increase of NT/PS, occurs. For 
larger values of e, the circulation is formed in a wider range of 0, as seen in Fig.[^. 

Since the circulation is formed above a certain threshold value of e, normalized total 
angular momentum M of the system is defined as an order parameter of the system. 

1 ^ 

M=—^T.i^^-Ro)Av.. (3) 

Here and are coordinates and velocities of i-th particle, Rq are the coordinates of the 
center of the box, "A" refers to outer product, and the notation v is used to remind us that 
we are observing the "R-system" . M is normalized by the box length to eliminate the 
dependency on the box size. 

Time developments of M for some values of e are shown in Fig.|^, where A^ = 256 and 
= 0.5. For e = 0, the value of M fluctuates around M = 0. As e is increased, the 
fluctuations of M increase. For e = 0.036, 0.038, the circulating mean flow is formed. The 



direction of the circulation sometimes turns over. As e increases further, the events of 
turning over become rare. In order to investigate the formation of the circulation in detail, 
the distributions of M are shown in Fig.Q because these distributions are symmetric around 
M = 0, only the region M > is shown. 

As denoted by Fig.^, it is found that the peak position of M, Mpeak becomes nonzero 
for e above a certain threshold value and increases continuously from zero with increasing 
e. We also observe a similar behavior when we vary the value of while fixing the value of 
e. Thus the circulation appears continuously for any direction in the parameter space(e, 0). 
In Fig.|, the distributions of M are broad in shape, which indicates that the value of M 
fluctuates around the value at the peak Mpeak- If the formation of the circulation is a phase 
transition accompanied by symmetry breaking in the thermodynamic limit (A^ — > oo), the 
width of the distribution will converge to zero in this limit. In order to confirm this scenario, 
the distributions of M are examined in the dependence on A^ in Fig.|^, where A^ is varied 
while A^e is kept constant. 

In Fig.^(a), it is found that Mpeak will converge in the limit, A^ oo. This convergence 
shows that Mpeak is a function of A^e for sufficiently large A^, because the value of A^e is fixed 
in Fig.^(a). 

From Fig.^b), it is found that the width of the fiuctuation of M decreases in the manner 
of 1/ \fN as A^ increases. Thus, we conclude that the circulation appears as a phase transition 
with symmetry breaking. The turning over events of the circulation observed in Fig.|^ are 
finite size effects. 

In Fig.^ Mpeak around the critical point is shown. The horizontal axis is A^e(l — 
where the second term 0{1/^/N) comes from finite size effects. This figure clearly shows 
that Mpeak is a function of A^e for sufficiently large A^. Further calculations are needed for 
the precise determination of the critical exponents. 

In this letter, we simulate inelastic disks in a square box and investigate the steady states 
of the system using a new observation frame wherein the energy of the system is virtually 
conserved by the rescaling of time. The important parameter of the system is A^e, where A^ 
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is the number of particles and e is inelasticity. At e = 0, the steady states are homogeneous 
equilibrium states. At sufficiently small A^e, the steady states are also homogeneous. For A^e 
above the critical points, the homogeneous steady states are no longer stable Then the 
phase transition with symmetry breaking occurs where the circulation appears continuously. 

The results shown in this letter are independent of the characteristics of the model. As 
noted before, the results are independent of the dispersity of the particle radii. The hard- 
core potential is not essential to the results. We have conffimed that the similar results are 
obtained in the case of soft-core potential: a simulation of the models using the discrete 



element method [T^ with Nose-Hoover thermostat [|TI| without the rescaling of time. We 
have also observed similar circulating mean flow in the system with inelastic hard spheres. 
Thus similar results would also be obtained in three dimensions. 

A number of problems remain to be investigated. For example, for larger e up to e = 1, 
are there one or more additional phases? Is hydrodynamic description by A^ ^ oo with 
fixing A^e really possible for any boundary conditions? Do fluid- and solid phases exist in 
e 7^ 0? How are the phenomena in the system investigated here related to those in the 
system with a real energy source, such as vibrating beds? 

It is noteworthy that the homogeneous states in the system with finite e are always 
unstable for sufficiently large A^, because the phase transition occurs at finite A^e. Thus the 
conservative system is a singular-limit system when we take the limit N ^ oo first. 

The author thanks S. Sasa and N. Nakagawa for their useful input and careful reading 
of the manuscript, N. Ito and Y.-h. Taguchi for their helpful advice, and Dept. of Math. 
Sci. at Ibaraki Univ. for their hospitality. The author also acknowledges the support from 
JSPS. 
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FIGURES 

FIG. 1. Dependency of pressure P on 0. The vertical axis is NT / PS and the horizontal axis is 
area fraction (j). (a) N = 256 and e = 0, 0.02, 0.04, 0.08. (b) N = 1024 and e = 0, 0.005, 0.01, 0.02. In 
the intermediate range of (j) and for sufficiently large e, differences of the dependency are remarkable. 
Comparing (a) with (b), a similar dependence is observed for the same values of A^e. 

FIG. 2. (a) Snapshot of the system for e = 0. (b) Snapshot of that for e = 0.08. Here, 
N = 256, = 0.5. The lines from the center of each particle show velocities of particles. Velocity 
correlation and a circulating mean flow are clearly observed in (b). 

FIG. 3. Time dependences of M. Here, N = 256, (j) = 0.5. The horizontal axes are time in the 
"R-system" (0 < i < 40000), while the vertical axes are M (-0.5 < M < 0.5). In each figure, the 
value of e is shown at the right side. For the parameters shown here, the mean free time is a value 
from 3.4 to 3.6. 

FIG. 4. Distribution of M. The horizontal axes are M (0 < Af < 0.5), where the vertical axes 
are frequency. The value of e is shown at the left hand side of each figure. Other parameters are 
same as in Fig.^. 

FIG. 5. A^-dependence of the distribution of M. Here, </> = 0.5 and A^e = 20.48. The number of 
particles are 64(dotted line) ,256(single dotted line),1024(broken line), and 4096(solid line), (a) 
The horizontal axis is M and the vertical axis is frequency, (b) The horizontal axis is the difference 
from the peak position times ^/N . The vertical axis is frequency which is normalized at the peak 
value. 

FIG. 6. Mpoak around the critical point. The vertical axis is Mpe^k ^ind the horizontal axis is 
r = Ne{l - 3/\/]V). The guide line is |r - 6.85|^/V6-5. 
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